Testing Navier-Stokes in package ViscousFlow

In [8]:

    

using ViscousFlow


    

In [9]:

    

using Plots
pyplot()
clibrary(:colorbrewer)
default(grid = false)


    

In [10]:

    

using Random


    

In [11]:

    

Random.seed!(1);


    

In [12]:

    

nx = 129; ny = 129;
w = Nodes(Dual,(nx,ny))
w .= rand(Float64,size(w));
Qq = Edges(Dual,w);
Ww = deepcopy(Qq);
ψ = deepcopy(w);


    

In [13]:

    

L = plan_laplacian(w,with_inverse=true)


    

    
Out[13]:





Discrete Laplacian (and inverse) on a (nx = 129, ny = 129) grid with spacing 1.0

In [15]:

    

@time nl = divergence(cellshift!(Qq,curl(L\w))∘cellshift!(Ww,w));


    

    




  0.095555 seconds (143.07 k allocations: 8.300 MiB)

In [16]:

    

woseen(x::Tuple{Real,Real},t;Re=1.0,x0::Tuple{Real,Real}=(0,0),t0=1) = 
                            exp(-((x[1]-x0[1])^2+(x[2]-x0[2])^2)/(4(t+t0)/Re))/(1+t/t0)


    

    
Out[16]:





woseen (generic function with 1 method)

In [17]:

    

Re = 200 + 50rand()
U = 1.0 + 0.2randn()
U∞ = (U,0.0)
Δx = 0.015;
Δt = min(0.5*Δx,0.5*Δx^2*Re);
xlim = (0.0,3.0);
ylim = (0.0,2.0);


    

In [18]:

    

sys = Systems.NavierStokes(Re,Δx,xlim,ylim,Δt,U∞ = U∞)


    

    
Out[18]:





Navier-Stokes system on a grid of size 202 x 136

In [19]:

    

w₀ = Nodes(Dual,size(sys));


    

In [20]:

    

xg,yg = coordinates(w₀,dx=Δx,I0=Systems.origin(sys))
x0 = (1.0,1.0); t0 = 1;
wexact(t) = [woseen((x,y),t;Re=Re,x0=x0.+U∞.*t,t0=t0) for x in xg, y in yg]


    

    
Out[20]:





wexact (generic function with 1 method)

In [21]:

    

ifrk = IFRK(w₀,sys.Δt,
                (t,w) -> Systems.plan_intfact(t,w,sys),
                (w,t) -> Systems.r₁(w,t,sys) ,rk=TimeMarching.RK31)


    

    
Out[21]:





Order-3 IF-RK integator with
   State of type Nodes{Dual,202,136}
   Time step size 0.0075

In [22]:

    

t = 0.0
w₀ .= wexact(t)
w = deepcopy(w₀)
tf = 1.0
T = 0:Δt:tf;


    

In [23]:

    

t = 0.0;
for ti in T
    global t, w = ifrk(t,w)
end


    

In [24]:

    

using LinearAlgebra


    

In [25]:

    

LinearAlgebra.norm(w-wexact(t),Inf)


    

    
Out[25]:





0.0011706795352210975

In [26]:

    

plot(xg,yg,w)


    

    
Out[26]:

In [27]:

    

plot(xg,w[:,65],label="numerical",ylim=(0,1))
plot!(xg,wexact(t)[:,65],label="exact")


    

    
Out[27]:

In [ ]: